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Abatraat  -  In  this  ~tvu  pei^  paper,  a  novel  procedure 
for  generating  an  ABMA  spectral  model  of  a  wide  sense 
stationary  time  series  Is  developed.  The  paraawters 
of  this  model  are  selected  so  chat  they  most  closely 
fit  a  sec  of  Tule-Walker  equations  which  are  estimated 
from  a  finite  sec  of  Cine  series'  observations.  This 
ARMA  modeling  method  has  been  found  Co  exhibit  a  spec¬ 
tral  estimation  performance  which  Is  typically  superior 
CO  such  altemaclves  as  Che  maximum  entropy  (AR)  method, 
classical  Fourier  procedures  (MA),and,  Che  Box- Jenkins 
method  (ABMA) . 

One  of  Che  principal  features  of  this  spectral 
estimation  method  Is  the  elegant  algebraic  structure  of 
the  linear  system  of  equations  which  need  be  solved 
when  finding  the  ARMA  model's  parameters.  This  shift- 
invariant  type  structure  gives  rise  to  an  adaptive 
algorithmic  solution  procedure  whose  computational 
efficiency  Is  comparable  to  that  achieved  by  recently 
developed  fast  AR  algorithmic  methods.  The  details  of 
Che  adaptive  ARMA  modeling  procedure  will  be  covered  In 
Part  2  of  this  paper.  These  dual  characteristics  of 
excellent  estimation  performance  and  real  time  adaptive 
Implementation  mark  this  method  as  being  a  primary 
spectral  estimation  cool. ^ 

I.  INTRODUCTION 

In  msny  Incerdlsclpllnary  applications,  it  Is  desired 
CO  estimate  the  essential  attributes  of  a  generally 
complex  valued  wlde-sense  stationary  time  series 
(x(n)).  Depending  on  the  specific  nature  of  Che  time 
series,  this  characterization  Is  often  adequately 
revealed  through  knowledge  of  the  time  series'  associa¬ 
ted  autocorrelation  sequence 

r^(n)  «  E{x(a+m)x*(m) )  n“0,  11.  i  2,...  (1) 

in  which  E  and  *  denote  the  operations  of  expectation 
and  complex  conjugation,  respectively.  On  Che  other 
hand,  the  requisite  characterization  may  often  be 
better  made  In  Che  frequency  domain  through  the  spect¬ 
ral  density  function 


r^(n)e 


which  Is  recognized  as  being  Che  Fourier  transform  of 
the  autocorrelation  sequence.  Either  member  of  this 
transform  pair  conveys  the  total  second-order  ststls- 
tlcal  Information  rslatlve  to  the  underlying  time 
series.  Frequently,  this  second  order  statistical 
characterization  provides  all  the  Information  required 
for  a  given  application  (e.g.,  optimal  Wiener  filtering, 
one-step  prediction,  etc.). 

The  classical  spectral  estimation  problem  la  con¬ 
cerned  with  estimating  the  spectral  density  function 
(2)  from  a  finite  sec  of  elms  sariss  observations. 
Without  loss  of  isnsrallty,  chess  observations  will  be 
taken  to  be  the  following  N  contiguous  elements 
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x(l),  x(2) . x(N)  s—  -  —  (3) 

A  variety  of  procedures  luve  been  proposed  for  using 
these  observations  to  effect  a  spectral  density  esti¬ 
mate.  Invariably,  Che  resultant  estimate  iflll  take  on 
a  rational  model  form  as  expressed  by 

bp,+b,e  .  .  .  +  b  e 

S  (eJ“)  -  -2 - i - 9 - 

*  1  +  a.e-^“+  .  .  .  ^  a  e’^P" 

1  P 

B  (e>)(^ 

-  9  V  (6) 

Ap(eJ“)| 

in  which  Che  ajc  and  bj^  are  referred  to  as  Che  model's 
autoregressive  and  moving  average  coefficients,  res¬ 
pectively.  We  shall  refer  to  this  particular  rational 
form  as  an  autoregressive-moving  average  (ARMA)  model 
of  order  (P>q).  It  Is  well  known  that  any  continuous 
spectral  density  can  be  approximated  arbitrarily 
closely  by  this  rational  form  if  the  order  pair  (p,q) 

Is  selected  adequately  large.  Thus,  by  Imposing  a 
rational  form  on  the  spectral  model,  we  Incur  no  real 
loss  In  spectral  representation. 

The  preponderance  of  research  and  application  inter¬ 
est  has  been  focused  on  two  special  cases  of  the  above 
ARMA  model.  They  are  the  moving  average  (MA)  model  in 
which  all  of  Che  aj^  coefficients  are  set  to  zero,  and, 
Che  autoregressive  (AR)  model  for  which  all  of  the  b^ 
coefficients  except  bg  are  set  to  zero.  The  spectral 
density  estimate  arising  from  a  HA  model  Is  seen  to 
possess  no  poles,  and  as  such  It  Is  frequently  referred 
CO  as  an  all-zero  model.  Using  simller  reasoning,  Che 
AR  model  Is  referred  to  as  an  all-pole  model,  and,  the 
general  ARMA  model  Is  referred  to  as  a  pole-zero  model. 

Classical  Fourier  approaches  [1]  and  the  perlodogram 
method  [2]  are  procedures  which  ultimately  provide  a 
MA  spectral  density  model.  Similarly,  the  maximum 
entropy  method  and  linear  predictive  coding  are  tech¬ 
niques  that  result  In  AR  spectral  density  models. 
Undoubtedly,  Che  primary  reasons  for  interest  in  speci¬ 
al  esse  MA  and  AR  models  lie  In  the  fact  chat  they: 

(1)  are  amenable  to  a  tractable  analysis,  (11)  typical¬ 
ly  provide  adequate  spectral  estimation  performance, 
and  (111)  give  rise  to  coefficient  selection  procedures 
which  ere  Implamsntable  by  computationally  efficient 
algorithms. 

Desplts  this  predisposition  towards  HA  and  AR  models, 
a  growing  Interest  In  ARMA  models  Is  evident  [3]-[9]. 
This  Is  In  recognition  of  the  fact  that  the  more 
general  ARMA  model  usually  provides  superior  spectral 
estimation  performance  while  at  the  same  time  require* 
fewer  model  parameters  to  achieve  that  behavior.  It  1* 
because  of  these  very  factors  chat  a  number  of  ARMA 
modeling  procedures  have  been  proposed.  These  include 
the  Box-Jsnklns  maximum  likelihood  method  [3),  whiten-  ^ 
Ing  filter  approaches  [4],  [5],  and,  more  recently, 
Cadzow's  high  performance  method  [6j-[9].  This  latter 
method  has  bean  found  to  provide  a  spectral  estimation 
performance  which  typically  excels  chat  obtained  from 
Its  MA,  AR,  and  ARMA  counterparts.  ^ 
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In  this  papsr,  ws  first  characterise  the  oodellng  of 
a  pure  ARHA  time  series.  An  analytical  procedure  Is 
presented  for  deteialnlng  the  underlying  aj^  and  b)^  co¬ 
efficients  In  which  the  tine  series'  actual  auto¬ 
correlation  eleoient  values  are  used.  This  Idealistic 
situation  then  provides  the  Justification  for  Intro¬ 
ducing  Che  high  performance  method  In  which  Che  ABMA 
model's  coefficients  are  estimated  from  time  series 
observations  and  not  from  autocorrelation  values.  It 
la  shown  Chat  Che  p  autoregressive  ai^  coefficients  are 
obtained  by  solving  a  consistent  system  of  p  linear 
equations.  Uhen  using  this  direct  approach,  Che  com¬ 
plete  set  of  time  series  observations  (3)  are  Incorpor¬ 
ated  to  effect  a  single  spectral  estimate  In  one 
computational  effort.  This  approach  Is  typically 
referred  to  as  "block  processing".  Horeover,  by  using 
the  generalized  Levinson  algorithm  [10]-[11],  It  Is 
possible  to  solve  the  above  mentioned  system  of  linear 
equations  in  a  computationally  efficient  manner. 

In  Part  2  of  this  paper,  a  recursive  procedure  Is 
developed  In  which  the  ARMA  model's  coefficients  are 
updated  as  each  new  time  series  observation  becomes 
available.  In  this  "tlme-update  processing"  mode,  an 
adaptive  form  of  spectral  estimation  Is  thereby 
achieved.  One  of- the  particularly  attractive  features 
of  this  time-updating  mode  Is  Its  computational 
efficiency.  Specifically,  the  p  autoregressive  co¬ 
efficients  (In  actuality  prediction  errors)  are 
optimally  updated  with  each  new  time  series  observation. 
The  number  of  multiplication  and  addition  coeputations 
required  In  this  updating  Is  of  the  order  p.  Thus, 

Che  computational  complexity  of  the  high  performance 
AKMA  method  Is  competitive  with  recently  developed 
"fast"  AS  methods,  but,  Its  spectral  estimation  per¬ 
formance  Is  typically  far  superior.  The  time-up^te 
node  la  particularly  attractive  In  those  situations  In 
which  Che  time  series  being  characterized  Is  a  long 
ongoing  process  and  one  wishes  to  generate  a  time 
evolving  sequence  of  spectral  estimates  In  a  real  time 
setting. 


II.  ABMA  TIME  SERIES:  PERFECT  MODELING 

In  this  section,  the  second-order  statistical  char¬ 
acter  Izatlon  of  an  ABMA  time  series  will  be  presented. 
This  characterization  will  play  a  central  role  in  the 
high  performance  spectral  estimation  procedure  chat  Is 
CO  be  developed  In  subsequent  sections.  The  time 
series  {x(n)}  Is  said  to  be  an  ABMA  time  series  of 
order  (p,q)  If  it  Is  generated  according  to  the  causal 
linear  recursive  relationship 

q  P 

x(n)  -  \  b.w(n-k)  -  \  a,x(n-k)  (5) 

k"0  k“I 

In  which  (w(n))  Is  a  zero  mean  white  noise  excitation 
whose  Individual  elements  have  variance  one.  It  is 
readily  shown  that  the  spectral  density  corresponding 
to  the  response  time  series  (x(n)}  is  given  by  ex¬ 
pression  (A).  Thus,  there  Is  seen  to  be  an  equivalen¬ 
ce  between  a  rational  spectral  density  model  and  the 
response  of  a  causal  linear  system  to  a  white  noise 
ejccltatlon. 

He  will  now  direct  our  attention  to  developing  a 
syacaamtlc  procedure  for  Identifying  the  recursive 
system's  autoregressive  coefficients  (l.e.,  the  Sk) 
and  moving  average  coefficients  (l.e.,  the  b)c)  from 
the  response  time  aeries'  autocorrelation  elements. 

It  will  be  beneficial  to  consider  separately  the  tasks 
of  Identifying  these  two  different  sets  of  coefficients. 

Autorsxresslve  Coefficient  Identification 

The  autoregressive  coefficients  can  be  determined 
directly  upon  examining  the  autocorrelation  character¬ 
ization  of  recursive  system  (5).  This  Is  achieved  by 


first  multiplying  both  sides  of  this  recursive  ex¬ 
pression  by  x*(n-m)  and  then  taking  the  expected  value. 
This  Is  found  to  result  In  the  well  known  Yule-Ualksr 
equations 

P 

[  a^T^(m-k)  -  for  “  2 

k-1 

where  It  Is  important  to  note  that  the  lag  parameter  a 
Is  hare  restricted  to  exceed  the  nuawrator  order  para¬ 
meter  q.  As  a  side  note,  the  Yule-Walker  equations 
will  Involve  the  moving  average  coefficients  b]^  In  s 
nonlinear  manner  for  lags  0  <  m  <  q.  The  characteristic 
equations  of  expression  (6)  provide  a  straightforward 
procedure  for  obtaining  the  ABMA  model's  a|^  auto¬ 
regressive  coefficients.  This  formally  entails  expres¬ 
sing  the  first  "t"  Yule-Halker  eqiiatlons  (l.e.,  q4-l  <  m 
<  qA-t)  In  the  following  matrix  format 


txCq) 

••  tjCq-p+i) 

r^(<rtl) 

r^(q+l) 

..  rjj(q-P+2) 

*2 

rx(q+2) 

■  - 

a 

P 

• 

rjj(q+t-l)  r^(q+t-2) 

..  rjj(q-p+t) 

. 

rjj(q+t) 

la  which  the  Integer  t  is  taken  to  be  equal  to  or 
larger  than  the  model's  denominator  order  (l.e.,  t  t  p). 
This  linear  system  of  equations  may  be  compactly 
expressed  as 


-4 


(8) 


where  R?_  Is  a  txp  autocorrelation  matrix,  rf  Is  a  t*! 
cp 

autoeorrelstloo  vector,  and,  Is  the  ARMA  model's  pxl 
autoregressive  coefficient  vector.  In  this  representation 
the  subscripts  t  and  p  are  appended  to  designate  the  num¬ 
ber  of  Yule-Halker  equations  being  used,  and,  the  ARMA 
model's  denominator  order,  respectively.  Similarly,  the 
superscript  q  depicts  the  ARMA  model's  numerator  order. 

To  obtain  the  ABMA  model's  autoregressive  co¬ 
efficients,  one  then  simply  solves  the  consistent  sys¬ 
tem  of  linear  equations  (8).  Valuable  Insight  relative 
to  rational  spectral  density  modeling  Is  provided  upon 
closer  examination  of  the  autocorrelation  matrix's 
(l.e.,  r4^)  algebraic  structure.  It  Is  convenient  to 

express  this  characterization  In  the  following  theorem. 

Theorem  I:  Let  {rx(k))  designate  the  auto¬ 
correlation  sequence  which  Is  associated  with 
an  ARMA  time  series  of  order  (p,q).  The 
corresponding  systam  of  t  linear  equations 
In  m  unknowns  as  specified  by 


R. 


(9) 


has  a  unique  solution  provided  that  mqp  and 
n>q  for  any  value  of  tap.  Moreover,  the  rank 
of  the  txa  matrix  R°_  is  given  by  min  (a,p,t) 

provided  that  n^q  and,  by  min  (m.t)  for  Oin<q. 

A  proof  of  this  theorem  will  not  be  given  here,  since 
these  results  are  Implicitly  doctaented  In  various 
textbooks  and  papers  dealing  with  time  series.  It  Is 
Important  to  note  that  even  If  one  has  perfect  auto¬ 
correlation  knowledge  of  an  ARMA  time  aeries,  the 
evaluation  of  the  associated  autoregressive  co¬ 
efficients  entails  a  determination  of  the  order  pair 
(p,q).  This  ordering  Information  is  Implicitly  con¬ 
tained  In  the  algebraic  structure  of  the  autocorrela¬ 
tion  matrix  R^  ,  and,  can  be  obtained  by  examining  this 
CB 
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struccur*  for  various  conblnaclona  of  che  nonnegscive 
Integers  n  and  a. 


B^(eJ“)B^*(e-*“')  -  Ap(e-*“’)C^*  (6-^“)  +  A^*  (e-^")C^ (e-*") 

(16) 


Moving  Average  Coefficient  Determination 

To  determine  the  coefficients  associated  with  the 
ABMA  time  series,  it  will  be  beneficial  to  Introduce 
the  causal  image  of  the  time  series'  autocorrelation 
sequence  as  defined  by 

r^^  (n)  •  r^(n)u(n)  --jr^COdCn)  (10) 

in  which  u(n)  and  S (n)  denote  the  standard  unit-step 
and  unlt-Rronccker  delta  sequences,  respectively.  The 
autocorrelation  sequence  may  be  recovered  from  its 
causal  image  by  ualng  che  complex  conjugate  symmetry 
property  of  autocorrelation  sequences  (l.e.,  r^(-n)  ■ 
T^(a)).  This  reconstruction  rule  cakes  the  form 

rj^(n)  •  (o)  + 

Upon  taking  che  Fourier  transform  of  relationship  (11), 
we  have  the  required  spectral  density  expression 

-  S^feJ-")  +  s;(eJ‘'-)* 

-  IRetS^"^  (e^“)l  (12) 

where  S^^(e'^'")  denotes  the  Fourier  transform  of  the 
causal  image  sequence  (r^'*'(n)}  ■ 

In  what  is  to  follow,  a  parametric  procedure  for 
representing  S^'*‘(eJ‘“)  (and  therefore  S^(eJ“))will  be 

given.  This  will  first  necessitate  che  Introduction 
of  the  auxiliary  saquenca 

P 

c(n)  •  r  "^(n)  +  I  a.  r  (n-k)  ,  0<nsmax(q,p)  (13) 
*  k“l  *  * 

in  which  the  causal  autocorrelation  elements  as  gene¬ 
rated  by  relationship  (10)  and  the  autoregressive  co¬ 
efficients  as  obtained  upon  solving  che  system  of 
equations  (8)  are  used.  According  to  Che  Yule-Walker 
aquations  (6)  and  Che  causal  image  definition  (10), 
it  is  seen  chat  this  auxiliary  sequence  is  identically 
zero  outside  the  indexing  range  0<n<max(q,p) .  With 
this  in  mind,  the  Fourier  transform  of  relationship 
(13)  is  naxt  taken  and  results  in 

C  (e>)  -  I  c(n)e-J“"  (14a) 

n-0 


-  (1  + 


P 

I 

n“l 


n  X 


•  A 

P 


(e'")s/(e>) 


(14b) 


in  which  a  ■  max(q,p).  Upon  solving  this  relationship 
for  S,('''(elw)  3iuj  substituting  this  solution  into  ex¬ 
pression  (12),  che  desired  ABMA  spectral  density  is 
obtained 


,  C  (eJ") 
S  (.J“)  -  -i-r- 

-  Ap(e>) 


C^*(eJ“) 


* 


A 


z 


(a^")C,(a>)  +  A^(.-l“)c/(.>) 

A  (eJ“)  A  *(aJ“) 

P  P 


(15) 


In  order  to  datermlna  the  ABMA  modal's  bf^  moving 
avarage  coafflclancs,  we  next  use  this  relationship  in 
conjunction  with  expression  (4)  to  obtain 


A  spectral  factorization  of  this  expression  will  then 
yield  che  prerequisite  b)(  coefficients  (assuming  a 
minimum  phase  BgCel*^)). 

In  summary,  the  spaecral  density  and  che  associated 
a|(  and  b^  coefficients  which  characterize  che  ABMA 
time  aeries  of  order  (P,q)  may  be  determined  by  follow¬ 
ing  the  four  step  procedure  as  outlined  in  Table  1. 

To  carry  out  this  model  identification  scheme,  it  is 
seen  chat  knowledge  of  che  order  pair  (P,q)  and  che 

q+p+1  autocorrelation  elements  r^fO) ,  rxd), _ >rx(d'*P) 

need  be  available. 


1.  Solve  relationship  (8)  for  Che  p  autoregressive 
a)(  coefficients.  This  will  require  setting  c>p. 


2.  Generate  che  auxiliary  sequence  c(n)  and  its 
Fourier  transform  using  expressions  (13)  and 
(14a),  respectively. 


3.  The  desired  spectral  density  is  then  given  by 
expression  (15). 


4.  Perform  a  spectral  factorization  of  che  poly¬ 
nomial  Bq(eJw)E^(aJw)  33  given  by  equation  (16) 
to  obtain  the  minimum  phase  choice  of  the  b)^ 
coefficients. 


Table  1:  Generation  of  che  spectral  density  and  the 
ABMA  model  parameters  associated  with  a  given 
sac  of  autocorrelation  values. 

III.  HIGH  PERFORMANCE  METHOD  OF  ABMA  SPECTRAL 
MODELING 

It  is  possible  to  adapt  many  of  cha  ideas  of 
Section  II  CO  achieve  an  ABMA  spectral  estimate  when 
only  cha  time  series  observations  (3)are  available 
(and  not  autocorrelation  values).  We  shall  again  treat 
separately  the  cases  of  autoregressive  and  moving 
average  coefficient  determination. 

Autoregressive  Coefficient  Estimation 

To  implement  the  autoregressive  coefflclanc  selacclon 
process  as  represented  by  relationship  (8)  it  will  be 
necessary  to  compute  appropriate  autocorrelation  esti¬ 
mates  from  che  given  sec  of  time  series’  obaervstlons. 
The  high  performance  ABMA  method  effects  Chase  esci- 
macca  in  the  guise  of  a  convenient  matrix  format  which 
lands  itself  to  a  particularly  efficient  computational 
raallzaclon  [6|-19].  In  particular,  the  autocorrela¬ 
tion  matrix  and  vector  required  in  expression  (8)  are 
estimated  according  to 

R^’p  -  Y^X  (17) 

r^’’  -  Y*x  (18) 

where  the  dagger  symbol  T  denotes  che  operation  of  com¬ 
plex  conjugate  transposition.  The  (N-p)xp  Toeplitz 


type 

matrix 

X  is  specified  by 

*x(p) 

x(p-l)  .  . 

.  x(l)  ■ 

x(p+l) 

x(p)  .  . 

.  x(2) 

X  - 

i 

x(N-l) 

x(N-2) 

x(N-p) 

while  Che  fN-p)«t  Toeplitz  type  matrix  Y  has  the  form 
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t 


x(p-q)  x(p-q-l)  .  .  .  x(p-q-t+l)' 
x(p-q+l)  x(p-q)  .  .  .  x(p-q-t+2) 


appropriate  autoregressive  coefficient  vector,  we  shall 
Introduce  the  following  quadratic  functional 


Y  - 


x(l»-q-l)  x{N-q-2)  .  .  .  x(H-q-t) 


and  X  Is  aji(N>p)xl  vector  given  by^ 
X-  (x(p+l),x(p<-2),  .  .  .  x(I»)]' 


(20)  *  ^  (25) 

In  which  A  Is  a  txt  posltlve-sealdeflnlte  diagonal  oatrlx 
with  diagonal  elenants  Chat  la  Introduced  In 
order  to  provide  one  with  the  option  of  weighting 
differently  the  various  error  vector  cooponencs.  It  Is 
a  slaple  aatter  to  show  that  an  autoregressive  co¬ 
efficient  vector  which  will  render  this  quadratic 

(21)  functional  a  nlniinuB  must  satisfy 


In  focnilatlng  matrix  Y,  we  have  used  Che  convention 
of  setting  CO  zero  any  elements  x(k)  for  which  k  lies 
outside  Che  observation  Index  range  I  l  k  l  N. 

If  the  autocorrelation  matrix  and  vector  estimates 
(17)  and  (18),  respectively,  are  subscltucad  Into  Che 
Yule-Walker  relationship  (8),  however.  It  la  generally 
found  that  Che  resultant  system  of  c  equations  In  Che 
p  autoregressive  coefficients  Is  Inconsistent  for  op. 
This  Is  due  to  Inevitable  Inaccuracies  in  the  auto¬ 
correlation  estimates,  and,  to  a  possible  Improper  ASHA 
model  order  choice.  In  any  case,  the  system  of 
equations  with  these  estimate  aubsclcuclona  will  give 
rise  Co  the  cxl  Yule-Walker  approximation  error  vector 
as  specified  by 

a  -  Y+X  a  +  Y^x  (22) 


Upon  taking  the  expected  value  of  e.  It  la  found  chat 
for  Che  ASHA  modeling  order  choice  In  which  q>p,  Chat 
this  expectation  results  In 


E{e(k))  •  (»-q-k) 


rx(q+k)  +  I  a^r^(q+k-m) 


ISkst 


(23) 


while  for  the  modeling  order  case  q<p  this  expectation 
produces 

'  r  p 


E{e(k)} 


(**-p)  I  «  r,(<I+k-tt) 

m“l 


(N-q-k) jr^(q+k)  +  I  a^r^(q+k-m) 


.  lSk<p-q 

,  p-q<k<t 
(24) 


In  either  ordering  case.  It  Is  seen  that  whan  the  clam 
series  Is  an  ARHA  process  of  order  (p,q),  Che  expected 
value  of  Che  error  vector  e  can  be  made  equal  to  zero 
by  a  proper  choice  of  the  autoregressive  coefficient 
vector  a.  Namely,  this  selection  would  be  such  chat 
the  underlying  Yule-Walker  equations  (6)  are  satisfied.^ 
This  Implies  chat  the  system  of  equations  (22)  with 
e  •  0  provides  an  unbiased  and  a  consistent  estimate  of 
the.  Yule-Walker  equations  (8) ,  where  £  is  the  zero  vector. 

With  the  above  thoughts  In  mind,  an  appealing 
approach  to  selecting  the  autoregressive  coefficient 
vector  Is  Immediately  suggested.  Namely,  a  Is  chosen 
so  as  to  make  the  arror  vector  "as  close"  to  its 
expected  value  of  ^  as  possible.  This  Is  of  course 
predlcsced  on  the  assumption  that  the  time  series  is  an 
AIKA  process  of  order  (p,q)  or  less.  In  order  to 
attain  a  tractable  procedure  for  selecting  an 


A  more  generalized  version  of  this  eatlawtlon  scheme 
can  be  obtained  by  substituting  the  Integer  k  for  p 
wherever  p  appears  In  relationship  (19)-(21).  For  ease 
of  presentation,  k  Is  here  restricted  to  be  p. 

2a  little  thought  will  convince  oneself  chat  this  same 
conclusion  will  be  reached  If  both  q  and  p  arc  at 
least  equal  to  the  niasarator  and  denoad.nator  orders, 
respectively,  of  the  underlying  ARMA  time  series. 


X+YAY'^Xa*  --X+YAY+x  (26) 

One  Chen  simply  solves  this  consistent  system  of  p 
linear  equations  In  Che  p  unknown  autoregressive  co¬ 
efficients  Co  obtain  an  estimate  for  Che  denominator 
of  Che  ASMA  model. 

Moving-Average  Coefficient  Estimation 

There  exist  several  procedures  for  estimating  Che 
ASMA  model’s  moving  average  coefficients.  We  shall 
now  briefly  describe  two  procedures  which  have  pro¬ 
vided  satisfactory  performance  and  In  a  sense  comple¬ 
ment  one  another. 


(1)  Cj^  Method 


The  procedure  which  has  provided  Che  best  fre¬ 
quency  resolution  behavior  Is  a  direct  adaption  of 
the  cv  method  as  described  in  Section  II  (see  ref. 
(8]).  In  particular,  using  Che  sec  of  auto¬ 
regressive  coefficient  estimates  as  obtained  from 
expression  (26)  and  a  suitable  sec  of  auto¬ 
correlation  estimates  fx(n)  for  nw0,l . max(q,p), 

one  computes  the  eje  coefficients  uai^  expression 
(13).  These  coefficients  are  then  used  to  achieve 
Che  desired  ASMA  spectral  estimate  when  incorporated 
Into  relationship  (14a)  and  ultimately  relationship 
(15).  Although  providing  an  excellent  frequency 
resolution  behavior,  this  procedure  suffers  the 
drawback  of  not  having  a  guaranteed  nonnegatlve 
definite  spectral  density  estimate^.  It  Is  with  this 
In  mind  chat  the  following  procedure  was  evolved. 

(11)  Smoothed  Periodogram  Method 

In  the  szuothed  periodogram  approach,  one  first 
computes  Che  so-called  "residual  time-series 
elements  according  to  the  relationship  (see  ref. [9]} 


P 

e(n)  “  x(n)  +  I  i.’x(o-k)  for  p<n<N  (27) 

k-1  * 


In  which  the  a|^*  autoregressive  coefficients  as  ob¬ 
tained  by  solving  expression  (26)  are  Incorporated. 
From  this  relationship,  It  Is  apparent  Chet  the 
following  spectral  density  expression  holds 


,  S  (e^*) 

S  (eJ“)  -  / 

*  lAVe^") 


(28) 


If  Sx(e^")  Is  to  correspond  to  an  ARMA  spectral  model 
of  order  (P>4,  It  Is  clear  that  a  order  HA 
spectral  estimate  for  the  residual  spectral  density 
S£(a-1")  must  be  obtained  and  Chen  substituted  Into 
relationship  (28).  Ths  smoothed  periodogram  has 
bean  found  to  be  a  useful  cool  for  this  purpose. 

In  ths  smooched  periodogram  method,  one  first 
partitions  Che  computed  residual  elements  (27)  Into 


^Thls  shorccoqU.ng  may  be  superficially  avoided  by 
caking  Che  absolute  value  of  Che  spectral  estlmaca. 


A 


L  segacncs  each  of  length  q+1  as  specified  by 

e.  (n)  «  e(ii+p+l+kd)  0<n<q  (29) 

Osk<I.-l 


where  "d"  Is  a  positive  Integer  which  specifies  the 
time  shift  between  adjacent  segments.  These  Indivi¬ 
dual  segments  will  overlap  If  d^q  and  will  perfectly 
partition  the  residual  sequence  when  d>q-*-l.  In  order 
to  Include  only  computed  elements,  the  relevant  para¬ 
meters  oust  be  selected  so  that  q+p+l+(L-l)dsN. 

Next  the  perlodogram  for  each  of  thee  L  segments  is 
taken  and  these  are  averaged  to  obtain  the  desired  q^^ 
order  smoothed  perlodogram,  that  la 

.  |2l 

-jun 


4  1  L-1 


k-0 


1 

q+1 


I  w(n)e,  (n)e 
n-0 


(30) 


J 


where  w(n)  is  a  window  sequence  that  is  normally 
selected  to  be  rectangular  (t.e.,  w(n)“l  for  0<n<q). 
The  required  ARMA  spectral  model  is  then  obtained  by 
substituting  this  approximation  into  relationship 
(28)  thereby  giving 


S  (eJ'-) 


3^.(6^““) 

:Ap*(eJ“)!2 


(31) 


It  is  readily  shotm  that  the  smoothed  perlodogram 
procedure  results  in  a  desired  nonnegative  qih  order 
MA  spectral  density  estimate.  Unfortunately,  its  fre¬ 
quency  resolution  capability  is  generally  not  of  Che 
sane  quality  as  chat  of  the  ck  method.^  On  the  ocher 
hand,  the  smooched  perlodogram  method  provides  more 
smoothly  behaved  spectral  estimates  which  contain  fewer 
spurious  effects. 

To  summarize,  Che  required  ARMA  spectral  model  is 
obtained  by  following  Che  systematic  procedure  out¬ 
lined  in  Table  2.  The  numerator  dynamic  estimation 
procedure  to  be  used  will  of  course  depend  on  Che  par¬ 
ticular  characteristic  being  sought  (e.g.,  frequency 
resolution,  smoothnese,  etc.). 


1.  Specify  values  for  the  ARMA  model's  order 

parameter  pair  (p.q),  the  Yule-Walker  equation 
parameter  c,  and,  Che  weighting  matrix's 
diagonal  elements 


2.  Using  Che  time  series  observations 

x(l) ,x(2) , . . . ,x(N) ,  construct  the  matrices 
X,  Y,  and  vector  x  according  to  relationships 
(19),  (20),  and  (21),  respectively. 

I - 

I  3.  Determine  the  model's  autoregressive  co- 
1  efficients  by  solving  relationship  (26) 


4.  The  numerator's  dynamics  are  obtained  by  using 
either  Che  (i)  c^  method,  or,  (11)  the  smooth¬ 
ed  perlodogram  method. 


Table  2.  Basic  steps  of  Che  standard  high  per¬ 
formance  ARMA  spectral  estimation  method: 

The  Block  Processing  Mode. 

The  improved  spectral  estimation  performance  ob¬ 
tained  in  using  this  high  performance  method  over  con¬ 
temporary  ARMA  techniques  such  as  Che  Box-Jenklns 
method  is,  to  a  large  extent,  a  consequence  of  select¬ 
ing  the  integer  c  to  be  larger  chan  Che  minimal 
value  p.  With  Che  corresponding  larger  set  of  Yule- 
Walker  equations  chat  are  thereby  being  approximated, 
it  intuitively  follows  chat  Che  model's  autoregressive 


coefficients  will  be  less  sensitive  Co  autocorrelation 
estimate  errors  which  are  embodied  in  Y'*3(  and  Y^x  Chan 
would  be  the  case  if  c  were  sec  to  p  (as  in  the  Box- 
Jenklns  method) .  This  anticipated  improvement  in 
spectral  estimation  behavior  when  using  Che  high  per¬ 
formance  method  has  In  fact  been  realized  on  a  rather 
large  number  of  numerical  examples  [6]-[9].  As  we  will 
see  In  part  2  this  high  performance  method  also  lends 
itself  to  a  particular  fast  adaptive  implemancaclon 
mode  when  C>p.  With  Che  two  attributes  of  Improved 
spectral  estimation  performance  and  coavutaclonal 
efficiency,  this  new  procedure  promises  to  be  an  import¬ 
ant  spectral  estimation  Cool. 

It  is  of  Interest  to  note  chat  when  q~0  and  t>p,  Che 
high  performance  ARMA  spectral  estimation  method  re¬ 
duces  Co  the  well  known  AR  covariance  method.  Moreover, 
upon  letting  c  exceed  p,  the  resultant  sec  of 
expanded  AR  Yule-Walker  equation  approximations  will 
typically  result  in  better  spectral  estimates  than  Che 
standard  AR  covariance  method.  To  Che  authors  know¬ 
ledge,  this  approach  has  not  been  used  in  Che  various  AR 
spectral  estimation  procedures  developed  to  dace. 


IV.  ORDER  SELECTION 

One  of  the  important  considerations  when  using  Che 
high  performance  method  is  chat  of  selecting  the  ARMA 
model  order  pair  (p.q).  This  selection  process  can  be 
made  by  utilizing  properties  of  the  ARMA  aucocorralac- 
ion  matrix  as  outlined  in  Theorem  1.  In  partlculsr, 
one  examines  the  column  rank  behavior  of  the  auto¬ 
correlation  matrix  estimate 


(32) 


chat  la  being  used  in  Che  high  performance  method. 

Upon  setting  q*t*p,  it  follows  that  the  p*p  auto¬ 
correlation  matrix  estimate  rP  will  start  becoming 
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ill-conditioned  when  the  order  parameter  p  exceeds  the 
time  series'  inherent  order  value  (assuming  that  q<p) . 
Thus,  the  model  order  determination  can  be  achieved  by 
investigating  the  conditioning  of  the  matrix  RP^  as  a 

function  of  p.  As  p  is  Increased,  an  appropriate 
choice  will  be  a  value  p  for  which  there  is  a  precipi¬ 
tate  decrease  in  matrix  conditioning  for  p-p  +  1.  This 
approach,  as  applied  to  Che  high  performance  method  of 
spectral  estlmaclon,  has  been  used  successfully  by  Pao 
and  Lee  [13] . 

There  exist  many  matrix  conditioning  measures  which 
may  be  used  for  this  order  determination.  One  of  the 
more  effective  measures  is  the  normalized  determinant 


as  specified  by  //p - - - ^ 

C(A)  -  det(A) //  [  I  la.,  r 
'  1-1  J-1 


(33) 


where  det(A)  designates  the  determinant  of  the  pxp 
matrix  A.  It  is  Co  be  noted  that  this  normalized 
determinant  will  be  zero  when  the  rank  of  A  Is  less 
chan  p. 


V.  THE  DOWN  SHIFT  OPERATOR 

In  the  analysis  to  follow,  extensive  use  of  the  dotm 
shift  operator  S  is  made.  This  operator  down¬ 

shifts  by  one  unit  the  elements  of  the  vector  upon 
which  it  operates  and  inserts  a  zero  into  the  vacated 
first  component  position.  In  ocher  words,  this  opera¬ 
tion  takas  the  form 

Sx  •  [0,  x(l),  x(2) . x(N-l)r  (34a) 

where  the  N'xl  vector  being  operated  upon  is  given  by 


^  t 

A  similar  approach  shares  Che  saisc  attributes  as  x  •  (x(l),  x(2) .  x(N)]  (34b) 

does  the  smooched  perlodogram.  [12] . 
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Th*  priai  syabol  hare  uaed  denoces  Che  operation  of 
vector  transposition.  It  Is  a  slsiple  aacter  to  show 
Chat  the  downshift  operator  has  Che  following  NxN 
natrix  representation 

S  ■  I«2  ;  ^3  ;  •  •  •  :  *1,  :  11  (35) 

in  which  ^  Is  Che  N^l  zero  vector  and  ^  dealgnaces  the 
kfh  standard  Nxl  basis  vector  %diose  components  are  all 
zero  except  for  Its  which  Is  one.  If  this  down¬ 
shift  operator  wore  applied  sequentially  o  times  to 
Che  vector  x.  It  is  clear  chat  a  downshift  of  o  units 
results,  chat  Is 

s“x  -  [0,  0 . 0,  x(l),x(2) . x(N-m)  (36) 


m  zeros 


VI.  PREWIMDOW  MODIFICATION 


In  many  spectral  estimation  applications.  It  is 
necessary  to  update  the  ARMA  model's  coefficients  as 
new  time  series  observations  become  available.  If  this 
Is  Co  be  achieved  in  real  time,  however.  It  Is  general¬ 
ly  not  feasible  to  apply  the  block  processing  Imple¬ 
mentation  of  the  high  performance  method  as  outlined  in 
Table  2.  In  Part  2  of  this  paper,  a  computationally 
efficient  algorithm  for  achieving  this  coefficient  up¬ 
dating  Is  developed.  In  order  to  facilitate  this  real 
time  recursive  algorithm.  It  Is  necessary  to  slightly 
modify  Che  constituent  matrices  X  and  Y,  and  Che 
vector  X  which  characterize  the  high  performance  method. 
These  modifications  provide  the  required  algebraic 
structure  to  render  the  resultant  modified  high  per¬ 
formance  ARMA  modeling  method  amenable  to  a  compu¬ 
tationally  efficient  recursive  solution. 

Although  a  number  of  modifications  are  possible,  we 
shall  only  treat  the  prewlndowlng  method  In  this 
Section. 3  in  the  premodlflcaclon  method,  the  x  vector 
Is  modified  to 

£-  [x(l),  x(2),  .  .  .,  x(N)l'  (37) 


while  the  X  matrix  is  modified  to  the  Nxp  Toepllcz 
type  matrix 


0 

x(l) 

x(2) 


0 

0 

x(l) 

x(2) 


x(N-l)  x(N-2) 


z(l) 


x(N-p) 


js*  :  s^x  .  .  .  :  s’’x  J 


(38) 


where  S  Is  the  downshift  operator.  Finally,  the  Y 
Mtrlx  Is  modified  Co  Che  Nxr  Toeplltz  type  matrix 


0 

0  ... 

0 

0 

0 

0 

0 

0 

x(l) 

0 

x(2) 

x(l) 

0 

x(l) 

i(M-q-l) 

x(N-q-2) 

x(N-q-t) 

^The  posewindowing,  and,  pre  &  postwindowing  modifica¬ 
tion  methods  are  described  In  Che  Appendix. 


(39) 


Upon  examination  of  these  expressions  for  Che  modified 
matrices  X  and  7,  It  Is  seen  chat  they  possess  a  very 
simple  shift  type  structure.  It  Is  this  very  structure 
which  renders  the  prewindowed  modification  amenable 
Co  a  coispuCatlonally  efficient  adaptive  solution  algor¬ 
ithm.  FurChenaore,  it  Is  to  be  noted  chat  lower 
triangular  pxp  and  pxt  matrices  have  been  added  to  Che 
cop  of  Che  original  X  and  Y  matrices  to  form  Che  modi¬ 
fied  X  and  7  matrices,  respectively.  These  augmenting 
lower  triangular  matrices  are  uniquely  specified  so  as 
to  make  the  modified  matrices  Toepllcz  In  structure 
(i.e.,  the  elements  along  any  diagonal  are  all  equal) 
with  zeros  appearing  In  the  upper  right  portion  of  each 
matrix.  It  is  this  specific  structure  wUch  makes  an 
efficient  recursive  solution  possible.  This  method  Is 
referred  Co  as  prewlndowlng  since  Che  Implicit  assump¬ 
tion  Chat  x(n)  >0  for  n<0  is  being  made. 

If  these  modifications  are  Incorporated  Into  expres¬ 
sion  (26),  a  modified  sec  of  p  linear  equations  In  Che 
p  autoregressive  coefficient  unknowns  Is  obtained,  chat 
Is 

X^y^Y*X  a  °  •  -X*itLY^x  (40) 

This  system  of  equations  represents  the  least-squares 
solution  CO  Che  following  statistical  approximation  of 
Che  first  c  Yule-Walker  equations 

e  -  7+7a  +7^  (41) 

The  effectiveness  of  this  approximation  can  be  evaluated 
by  caking  the  expected  value  of  this  relationship.  When 
the  ARMA  model  order  parameters  are  such  chat  q<p,  this 
expectation  Is  found  to  give 


E{e(k)) 


q+k  p 

(N-q-k)  I  a  r  (q+k-m)  +  J 

m"0  m^+k+1 


(N-m)a^rjj  (q+k-m), 
l<k<p-q 


P 

(N-q-k)  I  a  r  (q+k-m) 
m-0  "  * 


P*9Sk<t 

(42a) 


idiere  aQ  •  1.  This  Implies  chat  the  Yule-Walker  equat¬ 
ion  estimate  (42)  Is  biased  In  nature.  As  Che  data 
length  N  Increases,  however,  this  estimate  becomes 
asymptotically  unbiased.  For  Che  ordering  case  q>p,  the 
expectation  Is  found  to  yield 
P 

E{e(k)}»  (N-q-k)  7  a  r  (q+k-m)  ,  l<k<t  (42b) 

«H)  "  * 


which  is  unbiased  In  nature.  Thus,  Che  set  of  linear 
equation  estimates (41;  generally  provides  a  satisfactory 
estimate  for  the  associated  Yule-Walker  equations. 

In  order  to  achieve  Che  recursive  update  capability 
as  mentioned  previously,  it  will  be  necessary  to 
"restrict"  the  parasmter  t  to  be  p.  This  in  turn 
results  in  Z'^X  being  a  pxp  matrix.  When  this  matrix 
Is  Invertible,  there  always  exists  s  unique  auto¬ 
regressive  vector  %dilch  will  render  the  error  vector  to 
be  zero,  that  Is 

I*X  a*  -  -r^’x  (43) 

The  update  algorithm  to  be  presented  In  Part  2,  In 
effect,  allowa  us  to  recursively  obtain  Che  solution  for 
the  N+1  data  length  case  from  Che  solution  to  the  N 
data  length  case  [14] .  Unfortunately,  the  restriction 
of  c>p  also  generally  results  In  an  associated 
decrease  in  spectral  estimation  performance  (relative 
to  Op).  Thus,  In  obtaining  a  computationally  efficient 
update  recursive  algorithm,  an  accompanying  decrease  in 
spectral  estimation  performance  la  the  price  being  paid. 
One  must  therefore  carefully  consider  the  ramifications 
of  this  tradeoff  In  any  given  application.  It  Is  note¬ 
worthy,  however.  Chat  this  performance  degradation 
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diminishes  as  Che  number  of  time  series  observations  N 
grows. 

VII.  GENERALIZED  LEVINSON  ALGORITHM 

In  Che  high  performance  ASMA  modeling  procedures 
presented  In  Sections  III  and  VI,  Che  model's  p  auto¬ 
regressive  coefficients  were  obtained  by  solving  a  sys¬ 
tem  of  p  linear  equations.  In  the  special  case  In 
which  c>p  and  the  pxp  matrix  Y'*'X  la  nonslngular,  this 
relevant  system  of  equations  (26)  simplifies  to 

Y+Xa”  -  -Y+x  (44) 

where  Che  entries  of  Che  matrices  X  and  Y  and  the 
vector  X  are  dependent  on  Che  particular  form  being 
used  (l.e.,  unmodified,  prevlndowed,  poscwlndowed,  etc.) 

If  standard  matrix  Inversion  techniques  such  as  the 
Cholesky  decomposition  method  are  used,  on  the  order  of 
p3  mulcipllcaclona  and  additions  are  required  to  com¬ 
pute  Che  solution  to  relationship  (44).  These  standard 
techniques  ace  therefore  said  to  possess  a  computation¬ 
al  complexity  of  0(p3).  For  relatively  large  values 
of  p,  this  can  result  in  an  undesirable  computational 
burden.  On  the  other  hand.  If  the  pxp  matrix  has 
a  near  Toeplicz  structure.  It  is  possible  to  utilize 
the  generalized  Levinson  algorithm  to  obtain  the 
required  solution  using  far  fewer  computations  (10], 
[11).  Since  the  matrix  Y+X  Is  being  used  to  approxi¬ 
mate  the  Toeplitz  matrix  R*!  ,  there  is  good  reason  to 
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anticipate  that  Y  X  might  possess  this  structural 
feature. 

To  measure  Che  degree  to  which  Y^X  is  Toeplicz  In 
structure.  It  Is  necessary  to  introduce  the  concept  of 
displacement  rank.  The  displacement  rank  a (A)  of  Che 


pxp  matrix  A  Is  formally  given  by 

i(A)  -  min(a_(A),  a_|_(A)  ] 

(45a) 

where 

a_(A)  •  rank (A  -  SAS' ) 

(45b) 

a^(A)  •  rank(A  -  s'aS) 

(45c) 

In  which  S  Is  Che  aforementioned  down  shift 
(35).  When  the  matrix  A  is  Toeplitz,  It  is 

operator 

readily 

shown  that  Its  displacement  rank  is  two  (or  less). 

Thus,  a  matrix  whose  displacement  rank  Is  near  two  Is 
said  Co  be  close  to  Toeplitz  In  structure  and  therefore 
amenable  Co  efficient  inversion  using  the  generalized 
Levinson  algorithm. 

If  the  displacement  rank  of  Che  pxp  matrix  Y'X  Is  a. 
It  has  been  shown  chat  one  can  use  Che  generalized 
Levinson  algorithm  to  solve  expression  (44)  with  a 
corresponding  computational  complexity  of  0(up^)^.  If 
a  Is  sufficiently  smaller  chan  p,  a  significant  com¬ 
putational  savings  can  be  thereby  realized  relative  to 
standard  matrix  Inversion  routines.  Fortunately,  Che 
displacement  rank  of  Y'''X  Is  adequately  small  for  Che 
unmodified  high  performance  ARMA  modeling  method  and 
ICS  prewindowed  version  (as  well  as  the  poscwlndowed 
and  pre  4  poscwlndowed  versions).  This  Is  a  direct 
consequence  of  Che  fact  that  the  columns  of  matrices 
X  and  Y  are  simply  shifted  versions  on  one  another. 

One  may  readily  show  chat  the  dlsplacemanc  rank  of 
matrix  Y^X  for  each  of  the  high  performance  methods  Is 
as  shown  In  Table  3.  Sines  these  displacement  ranks 
are  so  small.  It  Is  clear  chat  the  generalized  Levinson 
algorithm  may  be  advantageously  used  for  solving  the 
linear  system  of  equations  (44). 


^As  a  byproduct  of  this  solution  procedure,  the  optimal 
autoregressive  coefficient  vectors  for  all  ARMA  models 
of  autoregressive  order  k  are  obtained  for  l<k<p. 


Method 

Displacement 

Rank 

a(YyX) 

Standard 

it 

Prevlndow 

3 

Postwindow 

3 

Pre  &  Postwindow 

2 

Table  3:  Displacement  rank  of  Che  matrix  Y'*'X 

for  the  various  high  performance  ARMA  methods 

When  Che  parameter  c  Is  allowed  to  Increase  beyond  p 
so  as  to  obtain  an  Improved  spectral  estimation  per¬ 
formance,  the  displacement  rank  of  each  of  the  methods 
spelled  out  in  Table  3  Increases.  It  is  readily  shown 
that  for  t>p  the  displacement  rank  increases  to 
a2(Y''’X)  in  all  cases.  For  example,  the  displacement 
rank  of  the  txp  matrix  Y'*’X  for  the  standard  procedure 
increases  to  (4)2  -le  and  so  forth.  For  excessively 
large  values  of  p,  it  would  then  be  advantageous  to  use 
the  generalized  Levinson  algorithm  to  solve  relationship 
(44)  when  case  t>p.  The  computational  complexity  there¬ 
by  obtained  would  be  on  the  order  of 

VIII.  numerical  example 

The  unmodified  ARMA  modeling  method  of  spectral  esti¬ 
mation,  as  presented  in  Section  III,  has  been  found  to 
possess  a  significantly  superior  performance  when  com¬ 
pared  to  such  contemporary  alternatives  as  the 
perlodogram,  maximum  entropy,  and,  the  Box-Jenklns 
methods  when  applied  to  "narrow"  band  time  series  (i.e., 
suosof  sinusoids  in  white  noise  [6l-[9]  and  [13]). 

With  this  in  mind,  the  effectiveness  of  both  the  un¬ 
modified  and  modified  ARMA  modeling  procedures  will  now 
be  examined  for  a  "moderstely  wide  band”  time  series. 

In  particular,  we  shall  treat  the  time  sarles  as  recent¬ 
ly  considered  by  Bruzzone  and  Ksveh  115).  Specifically, 
their  ARMA  time  series  of  order  (4,4)  is  characterized 
by  ,  7 

x^  *  ^  ^k  ^'^^k  (^^) 

1  2 

where  the  individual  time  series  x^  and  x^^  are  generated 
according  to 

\  ■  °-"Vl  -  °-”*k-2  ^  ^k 

x2. -0. 5x^-1  -0.93x^-2  +  t^ 

1  2 

In  which  the  Cy,  and  are  uncorrelated  Gaussian 

random  variables  with  zero  mean  and  unit  variance.  It 
then  follows  that  the  spectral  density  characterizing 
time  series  (46)  Is  given  by 

S  (w)  -  |l -0.4e''^“'+0.93s"-'^“i 

X  2 

•f  jl +0.5e'-)"  f0.93e"-l^“l  +0.25  (47) 

Using  the  time  series  description  (46),  twenty 
different  sampled  sequences  each  of  length  64  were 
generated.  These  twenty  observation  sets  were  then  used 
to  test  various  spectral  estimation  methods.  In  Figure 
1,  the  twenty  superimposed  plots  of  the  ARMA  model 
spectral  estimates  of  order  (4,4)  obtained  using  the 
first  Iterate  of  the  Box-Jenklns  method,  and,  this 
paper's  unmodified  method  with  •  (0.95)'^“3  and 
selections  of  t-4,  8,  and,  20  are  shown.  For  compari¬ 
son  purposes,  the  Ideal  spectrum  (47)  Is  also  shown. 

From  these  plots,  two  observations  may  be  made; 

(1)  the  unmodified  method  with  t-4  yields  a  marginally 
better  spectral  estimate  than  the  Box-Jenklns  method, 
and,  (11)  the  unmodified  spectral  estimates  Improve 
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significantly  as  c  Is  Increased  froa  the  nlnlaal 
value  4.  This  latter  observation  Is  most  noteworthy 
and  indicates  that  tha  Incorporation  of  sure  than  the 
■Inlnal  nuad>er  of  Yule-Walker  equations  for  detemln- 
Ing  the  AKMA  nodel's  autoregressive  coefficients  has 
the  anticipated  effect  of  significantly  iaprovlng 
spectral  estlaatlon  perfonsance. 

Next,  the  aodiflcatlon  aethods  developed  In  Section 
V  and  the  appendix  were  applied  to  these  twenty  differ¬ 
ent  saapled  sequences  of  length  64  to  obtain  ARMA  aodel 
spectral  estlaates  of  order  (4,4).  The  resultant 
spectra  are  shown  In  Figure  2  where  It  Is  apparent  that 
only  "a  modest"  degradation  In  spectral  estlaatlon  per- 
foraance  has  accrued  due  to  the  transient  effects  in¬ 
troduced  by  Che  modified  methods.  This  Is  indeed 
welcoaed  news  given  Che  ablllcy  co  lapleaent  chese 
modified  methods  with  exceptionally  fast  algorithas. 

It  Is  CO  be  noted  chat  the  "postaodlfled",  and  Che 
"pre  &  poscmodlfied"  methods  are  identical  in  this 
exaaq>le. 

As  a  final  example,  twenty  different  sampled  sequen¬ 
ces  each  of  length  200  were  generated  according  co 
expression  (46).  With  Chls  longer  data  length,  it  was 
anticipated  chat  an  improvement  in  spectral  estimation 
performance  would  result.  A  marked  improveaenc  is  in 
fact  realized  as  is  made  evident  from  Figure  3  where 
the  ARMA  model  spectral  estimates  of  order  (4,4)  are 
shown  for  the  Box-Jenklns  method  and  the  unmodified 
method  for  selections  of  t  ~4,  8,  and  20. 

IX.  CONCLUSION 

A  computationally  efficient  closed  fora  method  of 
ARMA  spectral  estimation  has  been  presented.  It  is 
predicated  on  the  approximation  of  a  sec  of  Yule-Walker 
equation  asclmaces  which  are  generated  froa  a  given  sec 
of  time  series  obsarvacloos.  The  ARMA  model's  auto¬ 
regressive  coefficients  are  determined  by  solving  a 
consistent  system  of  linear  equations.  The  displace¬ 
ment  rank  of  the  matrix  corresponding  to  these 
equations  is  four  thereby  indicating  that  an  efficient 
algorithmic  solution  procedure  is  possible. 

The  spectral  estimation  performance  of  this  ARMA 
modeling  procedure  has  been  eiaplrlcally  found  to  exceed 
chat  of  such  counterparts  as  the  ■««»<—■«  entropy  and 
Box-Jenklns  methods  (e.g.,  see  refs.  [6]-[9]  &  [13]). 
This  behavior  Is  Co  a  large  extent,  a  consequence  of 
Che  fact  that  more  chan  the  minimal  number  of  Yule- 
Walker  equation  estimates  are  being  approximated  Co 
obtain  the  resultant  ARMA  model  parameters. 

In  order  co  achieve  an  improved  computational 
efficiency,  a  prewindowed  modification  of  Che  proposed 
ARMA  model  spectral  method  was  next  Introduced.  The 
spectral  estimation  performance  of  this  prewlivdowed 
version  has  been  found  to  be  of  high  quality  for 
moderate  data  lengths.  As  ve  will  see  in  Part  2,  this 
prewindowed  method  may  be  laplenenced  by  an  adaptive 
update  algorlcha  whose  coaputaclonal  efficiency  Is 
comparable  to  chat  achieved  by  recently  developed  LMS 
fast  algorlchsia. 
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APPENDIX 


I.  Postwindow  Modification 


Following  a  similar  procedure  as  employed  In  Section 
VI,  the  addition  of  an  upper  triangular  matrix  to  the 
lower  portion  of  the  matrices  specified  by  equations 
(19)  aiul  (20)  yields  the  prevlndowed  matrices 


x(p)  .  . 

.  .  x(l) 

x(N)  . 

.  .  x(N-p+l) 

0  ■  ' 

•  x(N) 

(Al) 
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r 


x(p-q)  ....  x(p-q-t+l) 
x(p-q+l)  ....  x(p-q-t-«-2) 

z(N)  . x(N-t+l) 


O  x(H)  .  •  xINtp-q-t) 


(A2) 


where  X  and  7  are  recognized  as  being  (N^p)  and 
(Nxt)  Toeplltz  type  matrices,  respectively.  In  a 
similar  manner,  the  column  vector  x  is  nodlfied  to 

£"  [x(p+l) . .  x(N),  0,...,0]  (A3) 
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p  zeros 


Th«  displacamanc  rank  of  cht  matrix  T^X  la  raadlly 
found  CO  bo  y.  A  ganarallzad  Levinson  procedure  re¬ 
quiring  e  coapuCaclonei  complexity  of  OOp^)  can  then  be 
applied  for  solving  the  system  of  equations 

e  -  -y^z  {A4) 

A  more  computet lonelly  efficient  algorithm  assoclacad 
with  the  postvlnduw  modification  has  bean  developed  [141. 
It  la  shown  that  the  number  of  compucaclona  la  reduced 
to  (p  log  p)  If  p  •  q  ubara  p  and  q  are  the  denominator 
and  numerator  orders  of  the  ARMA  modal,  respectively. 


11.  Pre  &  Poscwlndou  Modification  Method 


The  combination  of  the  previously  discussed  pre- 
wlndovad  and  poscwlndowed  modification  methods  yields 
Che  pre  &  postvlndow  modification  method.  The  matrices 
and  vectors  are  modified  In  the  following  manner. 


u 

.  .  .  .  u 

x(l) 

■  •  0 

X  •  1  x(p) 

.  .  .■  x(l) 

! 

.  .  .  x(N-p) 

!  ^ 

u 

•  x(N) 

r?  •  • 

“l 

jo  .  . 

. 0 

i  X(l) 

.  ‘  ■  0 

i 

?  "  ' 

i 

x(t)  . 

.  .  .  .-.xd)  ; 

1  i(N)  . 

t+i)  j 

1  Q 

•x(H)  ...x(N+p-q-t)l 

C 

£  -  lx(l).  . 

..  x(N),  0,  . 

J 

..  01 

p  zeros 


(A5) 


(A6) 


(A7) 


where  X  and  X  denote  (M4-p)xp  and  {N+p)xt  Toeplltz  typo 
matrices,  respectively,  and  z  denotes  a  (Nfp)xl  column 
vector. 

tt  can  be  shown  chat  X^X  Is  a  Toeplltz  matrix.  The 
conventional  Levinson  algorithm  may  therefore  be  used 
for  solving  the  Toeplltz  system  of  equations 

i*;fa --y^z  (A8) 

in  which  Che  Inherent  computational  complexity  is 
0(2p-).^  More  recently,  a  fast  algorithmic  solution 
has  been  developed  which  significantly  reduces  this 
computational  complexity. 
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Figure  1;  Spectral  estimates  of  order  (4,4)  by 

the  Box-Jenklns  method  and  by  the  High 
Performance  method  using  various 
values  for  c.  N  ■  64  data  points  for 
each  estimate. 


The  parameter  t  is  here  taken  to  equal  p. 
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